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Abstract. In this paper we investigate the relationship between stabilized and enriched finite 
QQ ■ element formulations for the Stokes problem. We also present a new stabilized mixed formulation 

o 



for which the stability parameter is derived purely by the method of weighted residuals. This new 
formulation allows equal order interpolation for the velocity and pressure fields. Finally, we show by 
counterexample that a direct equivalence between subgrid-based stabilized finite element methods 
and Galerkin methods enriched by bubble functions cannot be constructed for quadrilateral and 
■ hexahedral elements using standard bubble functions. 
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1. INTRODUCTION 



£SJ . It is well known that under the mixed Galerkin formulation for the Stokes problem many prac- 
> 

Q^ ' tically convenient combinations of interpolation functions for the velocity and pressure fields often 

■ do not yield stable results. In particular, equal order interpolation for velocity and pressure (which 



is computationally the most convenient) is not stable. This numerical instability is attributed 
to the lack of stability in the pressure field which is mathematically explained by the celebrated 
Ladyzhenskaya-Babuska-Brezzi (LBB) stability condition pQ. One can verify the LBB condition 
numerically by means of a well designed patch test. 
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To address the deficiencies in the classical mixed formulation of the Stokes equations, two classes 
have emerged grouped by similar methodologies: stabilized finite element methods and enriched 
finite element methods. By stabilized finite element methods we mean methods that add mesh 
dependent terms to the standard Galerkin formulation that enable the formulation to satisfy or 
circumvent the LBB condition [2|. In contrast, enriched finite element methods add bubble functions 
to the finite element function space, which in turn play a stabilizing role. 
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Traditionally, the two most popular stabilized methods are the Streamline- Upwind/ 'Petrov-Galerkin 
(SUPG) method [3] and the Galerkin /least- squares (GLS) method [I]. In the GLS method, least- 
squares forms of the residuals are added to the Galerkin finite element formulation. These residual 
based terms are defined over the element interiors only, and the terms on the element boundaries 
are excluded. The underlying philosophy of the SUPG and the GLS methods is to strengthen 
the classical variational formulations so that the discrete approximations, which would otherwise 
be unstable, become stable and convergent. In the mid-90s Hughes revisited the origins of the 
stabilization schemes from a variational multiscale view point [5]. Under this variational multi- 
scale method, different stabilization techniques (including GLS and SUPG) are special cases of the 
underlying subgrid scale modeling concept. Although referred to as the Douglas- Wang method 
[6] rather than the variational multiscale method, Franca et al. [2] proved that the variational 
multiscale approach for Stokes flow is stable for many elements including T3, TET4, Q4, B8, and 
combinations of T3 and Q4 elements. 

Whereas the applicability of stabilized methods has been established for a wide range of elements, 
triangular elements have been the primary focus of enriched methods. Arnold et al. [7] were the 
first to develop the enriched finite element method for Stokes flow. Using continuous piecewise 
linear functions enriched by bubbles for velocity and piecewise linear functions for pressure, they 
showed that the MINI element satisfies the LBB condition. The stability of similar enrichment for 
quadrilateral elements has not yet been established. 

Although much work has been reported in the literature (see e.g., [8j (U \W\ lllj ) showing an 
equivalence between stabilized and enriched methods, the equivalence is not true for all PDE's nor 
is it true for all elements. For example, in the case of nearly incompressible elasticity the equivalence 
holds for T3, TET4, Q4, and B8 elements (see [12]). Likewise, for Stokes flow the equivalence holds 
for T3 and TET4 elements. However, we will demonstrate that the equivalence breaks down for 
Stokes flow for Q4 and B8 elements. 

The primary aim of this work is twofold: first, we present a new stabilized formulation for 
the Stokes problem, which can be derived purely by the method of weighted residuals. This new 
formulation is stable for equal order interpolation for the velocity and pressure fields, which is 
computationally convenient. Second, we show that enriching Q4 and B8 elements with standard 
bubble functions produces spurious pressure oscillations. The results confirm that one cannot 
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construct an equivalence between stabilized methods and bubble enrichment methods for these 
elements. 

The rest of this paper is organized as follows. First we present a consistent stabilized formulation 
for which the stability parameter is constructed from the element residual. Next, we propose 
an alternative residual based formulation that can be derived purely by the method of weighted 
residuals. Lastly, we present a mathematically equivalent enriched formulation and show that 
whereas the stabilized formulation does not show spurious pressure oscillations for a given test 
problem, the enriched formulation does. We conclude with some remarks regarding the equivalence 
between stabilized and enriched finite element methods for the Stokes problem. 

2. GOVERNING EQUATIONS FOR THE STOKES PROBLEM 

Let Q be a bounded open domain, and T be its boundary which is assumed to be piecewise 
smooth. Mathematically, T is defined as T := — Q, where is the closure of 0. Let the velocity 
vector field be denoted by v : $7 — > M. nd , where u nd" is the number of spatial dimensions. Let 
the (kinematic) pressure field be denoted by p : 0, — > R. As usual, T is divided into two parts, 
denoted by T v and r , such that n T* = and U T* = T. T v is the part of the boundary on 
which velocity is prescribed, and r is part of the boundary on which traction is prescribed. The 
governing equations for Stokes flow can be written as 

(1) -2vV 2 v + Vp = b in Q 

(2) V • v = in n 

(3) v = v p on T v 

(4) -pn + v(n ■ V)t> = t n on T* 

where V is the gradient operator, V 2 is the laplacian operator, b is the body force, v > is the 
kinematic viscosity, v p is the prescribed velocity vector field, t n is the prescribed traction, and n 
is the unit outward normal vector to T. Equation (pQ) represents the balance of linear momentum, 
and equation ^) represents the continuity equation for an incompressible continuum. Equations 
([3]) and dU) are the Dirichlet and Neumann boundary conditions, respectively. 

In the next section, we present the classical mixed formulation for the Stokes equations which 
will be the basis for the stabilized and enriched formulations. 
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3. CLASSICAL MIXED FORMULATION 

Before we present the classical mixed formulation for the Stokes equations, let us define the 
function spaces that will be used in the remainder of this paper. The function spaces for the velocity 
v(x) and the weighting function associated with velocity, denoted by w(x), are respectively defined 
as 

(5) V := {v | v G (H\n)) nd ,v = v p onT} 

(6) W := {w | w € (H 1 (p.)) nd , W! = 0on T v } 

where H 1 ^) is a standard Sobolev space pQ. In the classical mixed formulation the function space 
for the pressure p(x) and its corresponding weighting function q{x) are given by 

(7) V:={p\peL 2 (n)} 

where L 2 (Q) is the space of square-integrable functions on the domain Q. In the stabilized formu- 
lations, the function space for p(x) and q(x) will be defined as 

(8) V-={p\p^H 1 {ti)} 

For further details on function spaces refer to Brezzi and Fortin pQ. 

Remark 1. When Dirichlet boundary conditions are imposed everywhere on the boundary, that is 
r* = 0, the pressure can be determined only up to an arbitrary constant. In order to define the 
pressure field uniquely, it is common to prescribe the average value of pressure, 

(9) ( pdn = Po 

Jn 

where po is arbitrarily chosen (and can be zero). Then, the appropriate function spaces for the 
pressure that should be used instead ofV( defined in equation ([7j) ) is 

(10) Vo :={p\pe L 2 (Q), [ p dU = 0} 

Jn 

Another way to define the pressure uniquely is to prescribe the value of the pressure at a point, 
which is computationally the most convenient. 



The classical mixed formulation (which is based on the Galerkin principle) for the Stokes equa- 
tions can be written as: Find v{x) € V and p(x) £ V such that 

(11) a(w; v) + b(w; p) = f (w) VwGW 

(12) b(v;q)=0 Vq£V 

Let us define the bilinear forms as: 

(13) a(w;v) := Vw : 2uVv dU 

Jn 

(14) b(w;p):=- (V • w) p dft 

Jn 

and the linear functional as 

(15) f(w):= [ wbdU+ [ wt n dT 

Jn Jr* 

Once the weak formulation of the governing equations is established, the approximate solution 
based on the finite element method is determined in the usual manner. First one chooses the 
approximating finite element spaces, which (for a conforming formulation) will be finite dimensional 
subspaces of the underlying function spaces of the weak formulation. Let the finite element function 
spaces for the velocity, the weighting function associated with the velocity, and the pressure be 
denoted by V h C V, W h C W, and V h C V respectively. The finite element formulation of the 
classical mixed formulation reads: Find v h (x) £ V h and p h (x) £ V h such that 

(16) a{w h ;v h ) + b{w h ;p h ) = f(w h ) Vw h eW h 

(17) b{v h ;q h )=0 Vq h eV h 

For mixed formulations, the inclusions V h C V, W h C W, and V h C V are themselves not sufficient 
to produce stable results, and additional conditions must be met by these finite element spaces to 
obtain meaningful numerical results. A systematic study of these types of conditions on function 
spaces to obtain stable numerical results is the main theme of mixed finite elements. One of the 
main conditions to be met is the LBB inf-sup stability condition. 

Although the classical mixed formulation has many advantages (mainly its simplicity and ex- 
tensions to turbulent flows), it also has several numerical deficiencies. Most importantly, many 
combinations of shape functions for the velocity and pressure do not satisfy the LBB stability 
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condition and therefore exhibit unphysical oscillations in numerical simulations. As mentioned pre- 
viously, two classes of methods have been developed to overcome the limitations associated with 
the classical Galerkin approach; methods that augment the formulation with stabilizing terms to 
circumvent the LBB stability condition and those that enrich the function space to satisfy the LBB 
condition. 

4. VARIATIONAL MULTISCALE FRAMEWORK 

Hughes [5] proposed a variational framework based on the multiscale decomposition of the un- 
derlying fields into a coarse or resolvable scale and a subgrid or unresolvable scale. This framework 
provides a systematic procedure to develop stable finite element formulations. In this section, we 
present a multiscale formulation for the Stokes equations. A similar formulation for Darcy flow is 
presented in [13]. 

4.1. Multiscale decomposition. Let us divide the domain into N non-overlapping subdomains 
Q e (which in the finite element context will be elements) such that 

v 

(is) n=\Jn e 

e=l 

The boundary of the element O e is denoted by T e . We decompose the velocity field v(x) into 
coarse-scale and fine-scale components, indicated as v(x) and v'(x), respectively. To wit, 

(19) v(x) = v(x) + v'(x) 

Likewise, we decompose the weighting function w(x) into coarse-scale w(x) and fine-scale w'(x) 
components. 

(20) w(x) = w(x) + w'(x) 

We further make an assumption that the fine-scale components vanish along each element boundary. 

(21) v'(x) = w'(x) = on r e ; e = 1, . . . , N 

Let V be the function space for the coarse-scale component of the velocity v, and VV be the function 
space for w; and are defined as 

(22) V:= V; W := W 
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where V and W are denned earlier in equation ([5]) and equation ([6]) respectively. Let V' be the 
function space for both the fine-scale component of the velocity v' and its corresponding weighting 
function w', and is defined as 

(23) V :={v\v £ (H 1 ^))™ 1 , v = on T e , e = 1, . . . , N} 

The velocity field v{x) is now an element of the function space generated by the direct sum of 
V and V', denoted by V © V'. Similarly the direct sum of VV and V', denoted by VV © V', is the 
function space for the field w(x). 

In theory, we could decompose the pressure field into coarse-scale and fine-scale components. 
However, for simplicity we assume that there are no fine-scale terms for the pressure p(x) and for 
its corresponding weighting function q(x). Hence the function space for the fields p(x) and q(x) is 
V. 

4.2. Two-level classical mixed formulation. Substitution of equations (fl~9j) and (j20|) into the 
classical mixed formulation given by equations (jlip and (|12p becomes the first point of departure 
from the classical Galerkin formulation. 

(24) a(w + w'; v + v') + b(w + w'; p) = f(w + w') 

(25) b(v + v';q)=0 

Because the weighting functions w and w' are arbitrary, and because the functionals are linear 
in the weighting functions, we can write the above problem as two sub-problems. The coarse-scale 
problem can be written as: 

(26) a(w;v + v') + b(w;p) = f(w) V w € VV 

(27) b(v + v';q)=0 VgG? 

where the quantities a(-;-), &(•;■) and /(•;•) are defined in equations (fT3"j) - (fTSj) . The fine-scale 
problem can be written as: 

(28) a(w';v + v')+b(w';p) = f(w') Vw'eW' 

Remark 2. Note that the fine scale problem is independent and uncoupled at the element level 
(defined over the sum of element interiors). Due to the assumption that the subgrid scale response 
vanishes on the element boundaries, a(w;v + v') = a(w; v) + a(w; v'). 
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Using the linearity of the solution field and the divergence theorem on a(w;v'), we may alter- 
natively write the coarse-scale problem as: 

(29) a(w; v) + b(w;p) + c(w; v') = f(w) 

(30) b{v;q)+d(v';q) =0 

and the fine-scale problem as: 

(31) a(w'; v') + c(v; w') + d(w'; p) = f(w') 
where 

(32) c(w;v) := - / V 2 w • 2vv dQ 

Jn 

(33) d(w;p) := / w ■ Vp dQ 

5. FINE-SCALE INTERPOLATION AND BUBBLE FUNCTIONS 

If one chooses a single bubble function for interpolating the fine-scale variables (similar to the 
MINI element), then we have 

(34) v' = b e f3; w' = 6 e 7 

where b e is a bubble function, and f3 and 7 are constant vectors. The gradients of the fine-scale 
velocity and weighting functions are 

(35) W = (3Vb eT - Vw' = jVb eT 

where V6 e is a dim x 1 vector of the derivatives of the bubble function. Standard bubble functions 
for several elements are provided in Tabled! 

Table 1. Bubble Functions for Standard Finite Elements 



Element Bubble function 

T3 66(1-6-6) 

TET4 £166(1-6-6-6) 

Q4 (l_ e 2)(l-e 2 2 ) 

B8 (i_£)(i_g)(i_^) 



We shall substitute these expressions into the above subproblems in two different fashions, which 
brings us to the point of departure between stabilized and enriched methods. 



5.1. Weak variational multiscale formulation. In the spirit of a stabilized method, we elimi- 
nate the fine-scale variables by solving the fine-scale problem (equation (|3ip ) in terms of the coarse- 
scale variables. We then substitute the fine-scale solution into the coarse-scale problem (equation 
(|29p ) and solve the coarse-scale problem to obtain v(x) and p(x). This procedure also produces 
the familiar stabilization parameter, r, with which we augment the classical Galerkin formulation. 
Traditionally, one solves the fine-scale problem in terms of the coarse-scale variables in a weak or 
integral sense. For this reason, we refer to this method as the weak variational multiscale (WVM) 
formulation. 

5.2. Stabilization parameter. Typically, the stabilization parameter is derived in a consistent 
manner by incorporating the coarse-scale residual evaluated over the element. Examples of such 
formulations include the work of Masud and Khurram [14j for the Stokes equations and that of 
Nakshatrala et al [12] for nearly incompressible linear elasticity. The derivation proceeds as follows. 

Returning to equation (|31 1) , substituting equation (|34|) , and noting the arbitrariness of 7 we have 



where r := 2uV 2 v — Vp + b is the collection of the coarse-scale terms in the fine-scale problem. To 
solve for f3, one can make the approximation that in the limit of mesh refinement, the coarse-scale 
residual is constant over the element domain. Hence, r is moved outside of the integral in equation 
([36|) such that 



Remark 3. Note this is the only approximation introduced for this method, aside from the as- 
sumption that the subgrid scales vanish on the boundary (which is the key feature of the variational 
multiscale framework). 

Remark 4. In the case of T3 and TET4 elements, the statement that the coarse-scale residual is 
constant over the element domain is not an approximation, but is exactly true if b is constant. 



(36) 




(37) 
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Referring to equation (|34p . the fine-scale velocity may then be written as 

(38) v' = b e p = —rf 

2v 

where we have introduced the stabilization parameter r 



(39) 



IV6T dO 



-i 



5.3. Weak variational multiscale Galerkin formulation. Since we have an expression for the 
fine-scale velocity, we can substitute equation (|38|) back into the the coarse-scale problem to obtain 
a stabilized version of the Galerkin formulation 

(40) a(w;v) + b(w;p) - c(w; — r(2uV 2 v - Vp + b)) = f(w) VweW 

2v 

(41) b(v;q) - d(— t(2vV 2 v - Vp + b);q) =0 

2v 



Note the bilinear forms are defined in equations (|13 |) — (|15|) and (|32p . 

6. STRONG VARIATIONAL MULTISCALE FORMULATION 

We now present a new stabilized formulation for the Stokes problem that is consistently derived 
from the method of weighted residuals. Whereas traditionally the fine-scale problem is solved in 
a weak or integral sense, in the following formulation we solve the fine-scale problem in a strong 
sense. Therefore, we refer to this method as the strong variational multiscale (SVM) formulation. 

Using integration by parts and the linearity of the solution field, we may rewrite the fine-scale 
problem (given by equation pip ) as: 

(42) c(v'; w') + c(v; w') + b(w'; p) = f(w') 

Using the notation for the coarse-scale residual r := 2vV 2 v — Vp + b, the above equation can be 
written as, 



(43) / w' ■ (-2zA7V -r)dn = 



v 2 v' 

n 

Because w' is arbitrary and vanishes on the element boundaries and because v' is constrained to 
vanish on the element boundaries, the strong form of equation (|43|) is 

(44) 2uV 2 v' = -r in n e ;e = l...N 

(45) v' = on T e 
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Remark 5. The strong form may also be written as 

(46) C [v'] = -f(sc) in fl e ; v'{x) = on T e ; e = 1 . . . N 

where C [•] = 2vV 2 {-) is £/ie linear differential operator of the fine-scale problem. The analytical 
solution to equation (|46p over the element domain may be written as 

(47) v'(x) = - f G(x,y)f(y) dfl y 

where G(x, y) is the Green's function for the operator C. The potential for r to emanate from the 
element's Green's function has been pointed out in [5]. 

Obtaining an analytical solution for the Green's function that is valid for any element configu- 
ration is not always possible. Also, in order to get stable results an approximation to the Green's 
function will suffice. To this end, we approximate the solution using a single bubble function 

(48) v' = b e f3; V V = (vV)/3 

where V 2 6 e is defined as the Laplacian of the bubble function, which will never be zero (see 
Appendix). Substituting equation (j4~8l) into equation (j4"l"l) we have 

(49) (3 



2uV 2 b e 

We now have an expression for the fine-scale velocity v' 

b e 1 

(50) v = — =z r = Tf 

where r is the stabilization parameter defined as 

b e 

< 51 > -=v* 

A straightforward analysis shows that for an element with characteristic dimension h, the stabi- 
lization parameter r scales as h 2 . 

Remark 6. It is well-known in the mixed finite element literature (for example, see [UITSjj that r 
must scale as h 2 to guarantee convergence, which appears to be satisfied by (I5ip . 



Remark 7. The above stabilization parameter makes no approximations relative to the coarse-scale 

residual, as in the weak variational multiscale formulation. Therefore in the case of quadrilateral 

or hexahedral elements, no additional approximations are introduced preserving a mathematically 

exact correspondence with the enriched formulation presented below. 
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6.1. Weak variational multiscale Galerkin formulation. After substitution of equation (|50p 
into the coarse-scale problem (equations (|29p and (|30|) ) the resulting weak form is again expressed 
exactly as equations (|4T)j) and (jHJ). 

(52) a(w;v) + b(w;p) + c(w; ^-t(2vV 2 v - Vp + b)) = f(w) V w £ W 

(53) b{v;q) +cI(^t(2u\/ 2 v - Vp + b);q) =0 Vg££ 

Note that for linear elements like the T3 and TET4, V 2 w and X7 2 v will be exactly zero. 

7. ENRICHED FORMULATION 

For the enriched formulation we treat the coarse and fine-scale problems (equations ()26 j) — (|28p ) 
as two residual equations of the variables and p. Instead of analytically solving for v' in 

terms of the coarse-scale variables (as in a stabilized formulation), we use static condensation to 
solve the problem in a two stage manner. The emphasis in this section is placed on the solution 
strategy since it represents the most relevant features of the enriched formulation. 

7.1. Scalar residual. The scalar residual equations may be written as 

(54) r c (v; v',p) := a(w;v + v') + b(w;p) — f(w) 

(55) r p (v;v') := b(v + v';q) 

(56) rf(v;v',p) := a(w'; v + v') + b(w';p) — f(w') 
where the subscripts 'c', 'p', and 'f stand for coarse, pressure, and fine. 

7.2. Vector residual. To preserve the mathematical analogue to the strong variational multiscale 
formulation, we again choose a single bubble function for interpolating the fine-scale variables such 
that equation (|34p holds. As usual, v and its weighting function w may be expressed in terms of 
the nodal values v and w as 

(57) v = % T N T ; w = -[u T N T 

where AT" is a row vector of shape functions for each node. Substituting equations and (|34p and 

(|57p into equations (1540 - (I56p and noting the arbitrariness of w and 7, we can construct vector 
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residuals, R, that are the sum contributions of the vector residuals at the element level given as 



(58) 
(59) 
(60) 



R e c (v;v',p) :-- 
R e p (v;v') :- 
RJ(v;v', P ) :-- 



-I (N I)b dQ 



2u I B T vec[Vv + Vv'] dft - / vec[G]p 

if J fig 

- / N T V ■ (v + v') dVL 

Vv + W] dn - I B' T vec[I]p dft - j (b e I)b dtl 



2v I B >T vec 



To write more compactly, we have made the substitutions 



(61) 



B = G /; G := J T DN T 
B' = g /; g := V x b e 



where DN represents a matrix of the first derivatives of the element shape functions, J the element 
jacobian matrix, vec[-] is an operation that represents a matrix with a vector, and is the Kronecker 
product [16] (see Appendix). 

7.3. Stiffness matrix. Moving all applied force terms in R to the right hand side, we can write 
equations (|58p ~ (|60p in matrix form as 



(62) 



K 



V 




' fc ' 


V 






v' 




.ff. 



where / represents the sum of the element contributions to the applied forces, defined as 
(63) f e c = j (AT I)b dft; f e p = 0; f} = f (b e I)b 

J fig J S^e 

The global stiffness matrix, K, before static condensation, has the form 



dft 



(64) 



K 



K cc 


K C p 


Kef 


Kp c 


Kpp 


K P f 




Kfp 


Kff 
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where the element contributions are computed as follows 



(65) 



K' 

cc 


= 2v 


K' 


= -/ 




J^fc 


= 2v 



B T B dQ; K e cp 



vec[G]AT dft; K% f = 2v I B T B' dQ 



N T vec[G] r dn- K e =0; K 



B' T B dn; K 



fp 



g N dVt; K 



N T g T dfl 



) f = 2v 



b' t b' dn 



Using block Gauss elimination on equation (|62p . the fine-scale components can be condensed 
from the stiffness matrix. The resulting matrix equation can be written as 



(66) 



K 



V 




fc 


_ p _ 




.fp. 



The global stiffness matrix has the form 
(67) K 



CC -K c p 



where we have augmented the coarse-scale components with the fine-scale components at the ele- 
ment level as follows 



(68) 
(69) 
(70) 
(71) 



K cc = K e cc -K e cf (K} f r 1 K% 



K cp - K< cp- K lf{ K ff) lK )p 



Kpc — Kpc ~ Kpf(Kff) 1 K e j c 



K 



^ ~ K pf( K ff) lK % 



Similarly, the applied force vector has been augmented at the element level as 

(72) T c = r c -Kt f (K} f )^r f ; T P = -K; f (Kj f )-'r f 

After solving for the coarse scale variables from equation (|66p . the fine-scale variables can be 
recovered with post processing if desired. 



8. NUMERICAL RESULTS 

In this section we contrast the performance of the enriched formulation with that of the weak 
and strong variational multiscale stabilized formulations for various test problems. 
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8.1. Constant velocity and pressure problem. The constant velocity and pressure test prob- 
lem represents an extremely simple physical state, yet even the most sophisticated formulation must 
capture it without oscillations. The solution to the constant velocity and pressure test problem is 
v = (10.0, 0,0);p = 10.0, which by inspection satisfies the governing equations (equations CO)-®). 
The boundary conditions are defined in Figure Q] for the two-dimensional case. 

8.1.1. TET4 elements. As already mentioned in the introduction, the stability of the enriched 
method has been proved for triangular elements, but for the sake of completeness, we show that 
TET4 elements also perform well for the constant velocity and pressure problem. The results are 
shown in Figure El 

Remark 8. As an aside, the authors would also like to point out that for a well- centered triangle 
(WCT) mesh (triangles with no interior angles greater than or equal to 90 degrees), even the stan- 
dard Galerkin formulation produces no oscillations for the constant velocity and pressure problem. 
The results are shown in Figure A proof for the stability of such meshes is yet to appear in the 
literature. 

8.1.2. Q4 elements. As pointed out in Remark HJ the statement that the coarse-scale residual is 
constant over the element domain in the limit of refinement is exactly true for T3 and TET4 ele- 
ments for a constant body force, but in the case of Q4 and B8 elements, this statement is only an 
approximation. Due to the introduction of this approximation, the enriched and stabilized formula- 
tions produce starkly contrasting results for Q4 elements when applied to the constant velocity and 
pressure problem. Neither the weak or strong variational multiscale formulation shows any oscilla- 
tions in the pressure or velocity, but as shown in Figure HI one can see that the enriched formulation 
shows severe pressure oscillations. Brezzi and Pitkaranta [17] proposed a stabilizing technique to 
remedy such spurious modes by circumventing the LBB condition. To do so, one augments the 
enriched formulation with an added stability term e(Vq;Vp) where e ~ 0{h 2 ). This resolves the 
"missing" K pp term in the stiffness matrix before static condensation. Performing this augmenta- 
tion indeed weakly stabilizes the constant velocity and pressure problem, but this artificial term 
is not mathematically consistent. A similar approach, the Pressure Stabilizing/Petrov-Galerkin 
(PSPG) method [15] circumvents the LBB condition, but preserves consistency by applying a per- 
turbed weight function to all terms in the momentum equation. Although the PSPG method avoids 
oscillations in the pressure, the stability parameter is usually defined in an ad hoc manner. 
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In [181 f!9l [20l [2T] , the authors present an eigenvalue problem associated with the discrete LBB 
condition. Analysis of the eigenvalue spectrum reveals certain oscillatory modes, for example the 
pure pressure modes for which the associated eigenvalues are zero. The pure pressure modes consist 
of the hydrostatic mode and the checkerboard mode. The hydrostatic mode can easily be removed 
by properly prescribing the pressure boundary conditions, but the checkerboard mode is related 
to linear dependence in the discretized system of equations. Using a unit square discretized with 
a grid of n x n enriched Q4 and T3 elements, we present the results from a similar eigenvalue 
analysis in Figure [6J The results show that bubble enrichment removes the checkerboard mode for 
the T3 (MINI) element, but that the checkerboard mode remains for the enriched Q4 element. The 
presence of the checkerboard mode for the enriched Q4 element is consistent with the results shown 
in Figured! 

8.1.3. B8 elements. Results similar to the two-dimensional case are obtained when extended to 
three-dimensions. In particular, the B8 element also shows non-physical oscillations for this test 
problem that increase with mesh refinement. Figure [5] shows the results of the three-dimensional 
test problem for a coarse mesh and Figure [9] shows the results for a refined mesh. Notice that 
no oscillations are present for the results obtained with the weak or strong variational multiscale 
formulations as presented in Figures [7] and El 

8.2. Strong variational multiscale formulation. To further verify the strong variational mul- 
tiscale formulation of the Stokes problem, we present the results for some test problems along with 
a convergence analysis. 

8.2.1. Lid-driven cavity. The first problem evaluated is the well-known Lid-driven cavity problem. 
A description of the domain, along with the boundary conditions is shown in Figure [TUJ Contours 
of the velocity and pressure are shown in Figure [TTJ The results are in good accordance with other 
published results as shown in Table [2j 



Table 2. Position of the main cavity vortex 





element type 


y-location 


Present simulation 


B8 


0.753 


Donea and Huerta[22] 


Q2Q1 


0.756 
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8.2.2. Body force driven cavity. Another problem evaluated is the body force driven cavity taken 
from [22]. The problem geometry is the same as the Lid-driven cavity except that a velocity 
v x = v y = 0.0 is prescribed on the boundary and a constant body force is applied to the entire 
domain. The prescribed constant body force is given as 

bi = (12 - 24y)x 4 + (-24 + 48y)x 3 + (-48y + 72y 2 - 48y 3 + 12)x 2 

+ (-2 + 2Ay - 72y 2 + 48y 3 )x + 1 - Ay + 12y 2 - 8y 3 
b 2 = (8 - 48y + 48y 2 )x 3 + (-12 + 72y - 72y 2 )x 2 

(73) + (4 - 24y + 48y 2 - 48y 3 + 24y 4 )x - 12y 2 + 24y 3 - 12y 4 

The exact solution is 

v x = x 2 (l-x) 2 (2y-6y 2 +4y 3 ) 
Vy = - y 2 (\ - y) 2 (2x - 6x 2 + 4x 3 ) 

(74) p = x(l — x) 

Numerical results are shown in Figure [T2"| and they correspond well with other published results. 
The convergence properties of the strong variational multiscale formulation are shown in Figure 
[T3l To measure the error in the velocity, the L 2 norm is used, whereas the i? 1 -semi norm is used 
to compute the error in the pressure. Notice that the convergence rates are as expected for the 
Stokes problem using linear elements |15| . 



9. CONCLUSIONS 

We have introduced a new stabilized formulation for the Stokes problem that is appropriate 
for equal order interpolation for the velocity and pressure fields. The new formulation produces 
a scalar stabilization parameter that is consistently derived purely by the method of weighted 
residuals. We have also shown that an equivalence between enriched finite element methods and 
stabilized methods for the Stokes problem does not exist for certain elements. In particular, we 
have shown that enriching Q4 and B8 elements with standard bubble functions produces unstable 
results. Clearly, this work highlights the need for more emphasis in the development of bubble 
function enriched methods and the exact nature of their relationship to stabilized formulations. 
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APPENDIX 

Notation and definitions. Consider an n x m matrix A and a p x q matrix B 



aii ••• ai. 



fl n I ... Oy; 



61,1 ... b 



1,5 



6„ 1 ... 6; 



The Kronecker product of these matrices is an np x mg matrix, and is defined as 



AQB :-- 



a\ : \B . . . a\,rnB 
a n \B . . . a n m B 



The vec[-] operator is defined as 



vec[A] := 



• The divergence of the jacobian matrix. Consider the jacobian matrix J := a matrix 
of the coordinates of the nodes of an element, x, and the first and second derivatives of the shape 
functions DN and D 2 N, such that x = Nx, [DN] nm = and [D 2 N] nms = J^g . Starting 
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with a simple identity, we can derive the divergence of the jacobian matrix as follows 



JJ 



9 



J J- 1 



dx k 

9 Jim j— i ~ 9J mU 



mk 



dx k mk T ,Jzm dx 



j, 



.\dJi 



im j — 1 



pi dx k 



pm 



dx k 



k 

dx k 

dJ^k 
dx k 

dx k 

dx k 

dx k dx k 
V ■ J 1 



Sik 

d 

dx k 






(8 ik ) 



-J 



-1 9 



ij- 1 



p 1 dx k " mk 



d (9N n 



pi 



dx k \d& 



J 



i 

mk 



' d£ s 9^ 8x k 



ink 



- 7~ 1 r D 2 N T^ 1 T^ 1 

<J p i x nl u ^ nms'J mk J sk 



J 1 x T D 2 NJ 1 J T 



To further clarify, for a Q4 element, x, DN, and D 2 N are defined as 



x 



x\ yi 

%3 V3 



DN :-- 



dNi 
96 



gJVi 
96 



L dti 96 



d 2 N 1 d 2 N-j 



D 2 N :-- 



d 2 Nj d 2 N 4 



9696 9696 



d 2 N 1 
9696 



d 2 Nj 
9696 



• The Laplacian of a bubble function. Noting that /3 is a constant vector and making use of 
th divergence of the jacobian matrix as shown above, the Laplacian of a bubble function can be 

19 



computed as follows 



vV 



dxi dx^ ^ ^ 

d (d¥{i)di k 



03 



dxi V d£ k dxi) 3 



03 



dxi d^, k dxi ' d£ k dxi dxi J 3 

d db e (z) dt P dt k | db e {p d g& 

d(,p d£, k dxi dxi <9£ fc dxi dxi 



0j 



a db e (0 di P dz k db e (0 dik , d dN n au 36 



dip di k dxi dxi di k dx r d£ s d£ m dxi 
H be : J~ l J~ T - V,¥ T J~ 1 x T D 2 N : J~ l J~ T ) [5 



H b : J 1 J T - V i b eT J~ 1 x T D 2 N : J 1 J T 



03 



For trilinear B8 elements, D 2 N is denned as 



D N :-- 



d 2 N ± 9 2 iVi d 2 N 1 
9696 96 96 9^96 



d 2 N$ d 2 N$ d 2 N 8 
L 9696 9696 9696 



d 2 N 1 
9696 



d 2 N 8 
9696 



Note that D 2 N is a matrix representation of a third-order tensor. The matrix representing the 
second derivatives of the bubble functions, H b , is given as 



d 2 b e 



d 2 b e 



d 2 b e 



9696 9696 9696 



d 2 b e 
9696 
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Figure 1. Constant velocity and pressure test problem: domain and boundary conditions. 




(a) (b) 

Figure 2. Constant velocity and pressure test problem: (a) x- velocity and (b) 
pressure for 36 TET4 elements using the enriched formulation. 
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(a) (b) 

Figure 3. Constant velocity and pressure test problem: (a) x- velocity and (b) 
pressure for 336 well-centered triangle elements using a standard Galerkin formula- 
tion. 




X X 

(a) (b) 

Figure 4. Constant velocity and pressure test problem: (a) x- velocity and (b) 
pressure for 100 Q4 elements using the enriched formulation. 
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(a) 



(b) 



Figure 5. Constant velocity and pressure test problem: (a) x- velocity and (b) 
pressure for 8 B8 elements using the enriched formulation. 
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FIGURE 6. Eigenvalues, A, associated with the discrete Stokes problem for enriched 
Q4 vs. enriched T3 (MINI) elements for 1/h = 10 (a) close-up of pure pressure 
modes (b) all eigenvalues shown. 
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(a) (b) 

Figure 7. Constant velocity and pressure test problem: (a) x-velocity and (b) 
pressure for 8 B8 elements using the weak variational multiscale formulation. 




(a) (b) 

Figure 8. Constant velocity and pressure test problem: (a) x-velocity and (b) 
pressure for 8 B8 elements using the strong variational multiscale formulation. 
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Figure 9. Constant velocity and pressure test problem: (a) x- velocity and (b) 
pressure for 4096 B8 elements using the enriched formulation. 
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Figure 10. Lid-driven cavity: problem statement and boundary conditions. The 
non-leaky cavity approach is used here which resolves the discontinuity at the upper 
two corners of the domain by assuming that the corners belong to the vertical walls. 
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(a) (b) 

Figure 13. Convergence rates on a uniform mesh comparing the strong and weak 
variational multiscale formulations (a) L 2 norm of velocity (b) .f^-semi norm of 
pressure. 
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